16S WNV Antibiotics Analysis

Project Description: Examining the enteric bacterial microbiome of West Nile Virus infected mice (pre and post infection) to determine alterations in bacterial community strucutre or specific taxa which covary with lethality.

Primary Collaborator: Larissa Thakray (lthackray@wustl.edu)

Other relevant files: The following Phyloseq objects are available. Each is distinguished based on the 16S reference database used for taxonomic classification. RDP and Silva were processed through the species assignment workflow:

• ps0.wnv_antibiotics.gg.RDS (GreenGenes)
• ps0.wnv_antibiotics.rdp.RDS (RDP)
• ps0.wnv_antibiotics.silva.RDS (Silva)
• ps0.wnv_antibiotics.hitDB.RDS (hitDB)

References: * http://f1000research.com/articles/5-1492/v1

• http://benjjneb.github.io/dada2/tutorial.htm

Workflow details: The R commands below represent a full analysis of the following:

1. Sample quality control
2. ASV properties
3. Community Composition
4. Alpha diversity
5. Beta diversity

Load libraries

library("tidyverse")
packageVersion("tidyverse")

## [1] '1.2.1'

library("reshape2")
packageVersion("reshape2")

## [1] '1.4.3'

library("plyr")
packageVersion("plyr")

## [1] '1.8.4'

library("phyloseq")
packageVersion("phyloseq")

## [1] '1.23.1'

library("RColorBrewer")
packageVersion("RColorBrewer")

## [1] '1.1.2'

library("vegan")
packageVersion("vegan")

## [1] '2.4.5'

library("gridExtra")
packageVersion("gridExtra")

## [1] '2.3'

library("knitr")
packageVersion("knitr")

## [1] '1.18'

library("plotly")
packageVersion("plotly")

## [1] '4.7.1'

library("microbiome")
packageVersion("microbiome")

## [1] '1.0.0'

library("ggpubr")
packageVersion("ggpubr")

## [1] '0.1.6'

library("data.table")
packageVersion("data.table")

## [1] '1.10.4.3'

library("mgcv")
packageVersion("mgcv")

## [1] '1.8.22'

library("pairwiseAdonis") #https://github.com/pmartinezarbizu/pairwiseAdonis
packageVersion("pairwiseAdonis")

## [1] '0.0.1'

Read in data

# Load Phyloseq Object
# Selected RDP due to it's up-to-date nature and conservative taxonomy. Other files are also valid for anlysis but are not explored here
ps0 <- readRDS("./Data/ps0.wnv_antibiotics.rdp.RDS")
ps0

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 3229 taxa and 535 samples ]
## tax_table()   Taxonomy Table:    [ 3229 taxa by 7 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 3229 tips and 3227 internal nodes ]

# Load mapping file
map <- import_qiime_sample_data("./Data/mapping_wnv_antibiotics.txt")
dim(map)

## [1] 520  24

ps0 <- merge_phyloseq(ps0, map)
ps0

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 3229 taxa and 520 samples ]
## sample_data() Sample Data:       [ 520 samples by 24 sample variables ]
## tax_table()   Taxonomy Table:    [ 3229 taxa by 7 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 3229 tips and 3227 internal nodes ]

# Perform a few sanity checks
sample_variables(ps0)

##  [1] "X.SampleID"           "BarcodeSequence"      "LinkerPrimerSequence"
##  [4] "Name"                 "MiSeq"                "DNAplate"            
##  [7] "Well"                 "Sample.ID"            "Description"         
## [10] "X16SPlate"            "SampleNo"             "Run"                 
## [13] "DNAcell"              "Sex"                  "CageBeforeTreatment" 
## [16] "MouseBeforeTreatment" "CageDuringTreatment"  "DaysTreatment"       
## [19] "Day"                  "Treatment"            "Virus"               
## [22] "SurvivalStatus"       "DayDeadPostInfection" "Description.1"

ntaxa(ps0)

## [1] 3229

rank_names(ps0)

## [1] "Kingdom" "Phylum"  "Class"   "Order"   "Family"  "Genus"   "Species"

get_taxa_unique(ps0, "Phylum")

## [1] "Bacteroidetes"             "Proteobacteria"           
## [3] "Firmicutes"                "Verrucomicrobia"          
## [5] "Tenericutes"               NA                         
## [7] "Actinobacteria"            "Cyanobacteria/Chloroplast"
## [9] "Fusobacteria"

Factor reordering and renaming

# Remove Day -14 cohoused data
levels(sample_data(ps0)$DaysTreatment)

## [1] "D.14" "D0"   "D13"  "D16"  "D18"  "D20"  "D3"   "D7"

ps0 <- subset_samples(ps0, DaysTreatment != "D.14")
levels(sample_data(ps0)$DaysTreatment)

## [1] "D0"  "D13" "D16" "D18" "D20" "D3"  "D7"

# Remove uninfected samples
levels(sample_data(ps0)$Virus)

## [1] "Uninfected" "WNV2000"

ps0 <- subset_samples(ps0, Virus == "WNV2000")
levels(sample_data(ps0)$Virus)

## [1] "WNV2000"

# Remove taxa no longer part of the count table due to sample removal
summary(taxa_sums(ps0))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0       0       3    2353      10  928695

ps0 <- prune_taxa(taxa_sums(ps0) > 0, ps0)
summary(taxa_sums(ps0))

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##      2.0      4.0      9.0   4233.2     72.5 928695.0

# Factor re-ordering, relabelling, etc.
# Reorder Time points
levels(sample_data(ps0)$DaysTreatment)

## [1] "D0"  "D13" "D16" "D18" "D20" "D3"  "D7"

sample_data(ps0)$DaysTreatment <- factor(sample_data(ps0)$DaysTreatment, levels = c("D0", "D3", "D7", "D13", "D16", "D18", "D20"))
levels(sample_data(ps0)$DaysTreatment)

## [1] "D0"  "D3"  "D7"  "D13" "D16" "D18" "D20"

# Reorder Treatments
levels(sample_data(ps0)$Treatment)

## [1] "Amp"      "AmpMetro" "Metro"    "Vehicle"

sample_data(ps0)$Treatment <- factor(sample_data(ps0)$Treatment, levels = c("Vehicle","Metro","Amp","AmpMetro"))
levels(sample_data(ps0)$Treatment)

## [1] "Vehicle"  "Metro"    "Amp"      "AmpMetro"

# Relabel Treatments
sample_data(ps0)$Treatment <- factor(sample_data(ps0)$Treatment, labels = c("Vehicle","Metro","Amp","Amp + Metro"))
levels(sample_data(ps0)$Treatment)

## [1] "Vehicle"     "Metro"       "Amp"         "Amp + Metro"

ASV summary statistics

# Create a new data frame of the sorted row sums, a column of sorted values from 1 to the total number of individuals/counts for each ASV and a categorical variable stating these are all ASVs.
readsumsdf = data.frame(nreads = sort(taxa_sums(ps0), TRUE), 
                        sorted = 1:ntaxa(ps0),
                        type = "ASVs")

# Add a column of sample sums (total number of individuals per sample)
readsumsdf = rbind(readsumsdf,
                   data.frame(nreads = sort(sample_sums(ps0), TRUE),
                              sorted = 1:nsamples(ps0),
                              type = "Samples"))

# Make a data frame with a column for the read counts of each sample for histogram production
sample_sum_df <- data.frame(sum = sample_sums(ps0))

# Make plots
# Generates a bar plot with # of reads (y-axis) for each taxa. Sorted from most to least abundant
# Generates a second bar plot with # of reads (y-axis) per sample. Sorted from most to least
p.reads = ggplot(readsumsdf, aes(x = sorted, y = nreads)) +
  geom_bar(stat = "identity") +
  ggtitle("ASV Assessment") +
  scale_y_log10() +
  facet_wrap(~type, scales = "free") +
  ylab("# of Reads")

# Histogram of the number of Samples (y-axis) at various read depths
p.reads.hist <- ggplot(sample_sum_df, aes(x = sum)) + 
  geom_histogram(color = "black", fill = "firebrick3", binwidth = 150) +
  ggtitle("Distribution of sample sequencing depth") + 
  xlab("Read counts") +
  ylab("# of Samples")

# Final plot, side-by-side
grid.arrange(p.reads, p.reads.hist, ncol = 2)

# Basic summary statistics
summary(sample_sums(ps0))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##      86    2064   19679   24122   37592   97513

Sample assessment

# Format a data table to combine sample summary data with sample variable data
ss <- sample_sums(ps0)
sd <- as.data.frame(sample_data(ps0))
ss.df <- merge(sd, data.frame("ASVs" = ss), by ="row.names")

# Plot the data by the treatment variable
y = 1000 # Set a threshold for the minimum number of acceptable reads. Can start as a guess
x = "DaysTreatment" # Set the x-axis variable you want to examine
label = "Sample.ID" # This is the label you want to overlay on the points that are below threshold y. Should be something sample specific

p.ss.boxplot <- ggplot(ss.df, aes_string(x, y = "ASVs")) + 
  stat_boxplot(geom = "errorbar", position = position_dodge(width = 0.8)) +
  geom_boxplot(outlier.colour="NA", position = position_dodge(width = 0.8), alpha = 0.2) +
  scale_y_log10() +
  facet_wrap(~Treatment) +
  geom_hline(yintercept = y, lty = 2) + 
  geom_point(position=position_jitterdodge(dodge.width = 0.8), aes_string(color = "SurvivalStatus"), size = 1.2) +
  geom_text(data = ss.df, aes_string(x, y="ASVs", label=label), size=2) # This labels a subset that fall below threshold variable y and labels them with the label variable
p.ss.boxplot

write.table(ss.df, file = "./Results/asv_stats.txt", sep = "\t")

List of samples selected as outliers: c(“D16.M5”, “D16.M2”, “D18.M3”, “D18.M5”, “D18.K1”, “D20.M5”, “D20.M2”)

# Outlier samples: c("D16.M5", "D16.M2", "D18.M3", "D18.M5", "D18.K1", "D20.M5", "D20.M2")
nsamples(ps0)

## [1] 315

ps1 <- ps0 %>%
  subset_samples(
    Sample.ID != "D16.M5" &
    Sample.ID != "D16.M2" &
    Sample.ID != "D18.M3" &
    Sample.ID != "D18.M5" &
    Sample.ID != "D18.K1" &
    Sample.ID != "D20.M5" &
    Sample.ID != "D20.M2"
)
nsamples(ps1)

## [1] 309

saveRDS(ps1, file = "./Results/ps1.RDS")

ggpaired(ss.df, x = "DaysTreatment", y = "ASVs", outlier.shape = NA, id = "CageDuringTreatment", line.size = 0.5, line.color = "gray") +
  scale_y_log10() +
  facet_grid(~Treatment) +
  stat_compare_means(ref.group = "D0", hide.ns = TRUE, label = "p.signif") +
  theme(axis.text.x = element_text(size = 8)) +
  theme(axis.text.y = element_text(size = 8)) +
  theme(axis.title.x = element_blank()) +
  ylab("Read Counts")

##Overall sample relationship to evaluate sample outliers

Taxon cleaning

# Begin by removing sequences that were classified as either mitochondria or chlorplast
ntaxa(ps1) # Check the number of taxa prior to removal

## [1] 1795

ps1 <- ps1 %>%
  subset_taxa(
    Family  != "mitochondria" &
    Class   != "Chloroplast"
  )
ntaxa(ps1) # Confirm that the taxa were removed

## [1] 1372

Data transformations

Subsetting

Community composition plotting

# Create a data table for ggploting
ps1_phylum <- ps1 %>%
  tax_glom(taxrank = "Phylum") %>%                     # agglomerate at phylum level
  transform_sample_counts(function(x) {x/sum(x)} ) %>% # Transform to rel. abundance (or use ps0.ra)
  psmelt() %>%                                         # Melt to long format for easy ggploting
  filter(Abundance > 0.01)                             # Filter out low abundance taxa

# Convert Sample No to a factor because R is weird sometime
ps1_phylum$SampleNo <- as.factor(ps1_phylum$SampleNo)

# Plot - Phylum
p.ra.phylum <- ggplot(ps1_phylum, aes(x = SampleNo, y = Abundance, fill = Phylum)) + 
  geom_bar(stat = "identity", width = 1) +
  facet_wrap(Treatment~DaysTreatment, scales = "free_x", nrow = 4, ncol = 7) +
  theme(axis.text.x = element_blank()) +
  theme(axis.title.x = element_blank()) +
  ggtitle("Abundant Phylum (> 1%)")
p.ra.phylum

# Note: This is a nice place to output tables of data that you may want to use for other analysis, or to include as supplemental data for publication
# You can rerun the first bit of code in this chunk and change Phylum to Species for a table with all possible classifications
write.table(ps1_phylum, file = "./Results/phylum_relab.txt", sep = "\t")

ggplotly(p.ra.phylum)

# agglomerate taxa
glom <- tax_glom(ps1.ra, taxrank = 'Phylum')

# create dataframe from phyloseq object
dat <- as.tibble(psmelt(glom))

# convert Phylum to a character vector from a factor because R
# dat$Phylum <- as.character(dat$Phylum)

# Reorder Phylum levels from most -> least abundant
levels(dat$Phylum)

## [1] "Actinobacteria"  "Bacteroidetes"   "Firmicutes"      "Fusobacteria"   
## [5] "Proteobacteria"  "Tenericutes"     "Verrucomicrobia"

dat$Phylum <- factor(dat$Phylum, levels = c("Bacteroidetes", "Firmicutes", "Proteobacteria", "Tenericutes", "Actinobacteria", "Verrucomicrobia"))
levels(dat$Phylum)

## [1] "Bacteroidetes"   "Firmicutes"      "Proteobacteria"  "Tenericutes"    
## [5] "Actinobacteria"  "Verrucomicrobia"

levels(dat$Treatment)

## [1] "Vehicle"     "Metro"       "Amp"         "Amp + Metro"

dat$Treatment <- factor(dat$Treatment, levels = c("Vehicle", "Metro", "Amp", "Amp + Metro"))
levels(dat$Treatment)

## [1] "Vehicle"     "Metro"       "Amp"         "Amp + Metro"

p.boxplot.phylum.1 <- ggpaired(dat, x = "DaysTreatment", y = "Abundance", outlier.shape = NA, id = "CageDuringTreatment", line.size = 0.5, line.color = "lightgray", point.size = 0.3, width = 0.4) +
  #geom_point(size = 1) +
  #geom_jitter(width = 0.2, alpha = 0.7) +
  ylab("Relative Abundance") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  stat_compare_means(label = "p.signif", method = "t.test", ref.group = "D0", hide.ns = TRUE, label.y.npc = 1) +
  facet_grid(Treatment~Phylum) +
  theme(axis.text.x = element_text(size = 6)) +
  theme(axis.text.y = element_text(size = 6)) +
  theme(axis.title.y = element_text(size = 6)) +
  theme(axis.title.x = element_blank()) +
  theme(strip.text = element_text(size = 6)) +
  ylim(0,1.2)
p.boxplot.phylum.1

# Reduced to most abundant phylum
summary(dat$Phylum)

##   Bacteroidetes      Firmicutes  Proteobacteria     Tenericutes 
##             309             309             309             309 
##  Actinobacteria Verrucomicrobia            NA's 
##             309             309             309

dat.1 <- filter(dat, Phylum %in% c("Bacteroidetes", "Firmicutes", "Proteobacteria", "Tenericutes"))
dat.1 <- droplevels(dat.1)
summary(dat.1$Phylum)

##  Bacteroidetes     Firmicutes Proteobacteria    Tenericutes 
##            309            309            309            309

levels(dat.1$Treatment)

## [1] "Vehicle"     "Metro"       "Amp"         "Amp + Metro"

dat.1$Treatment <- factor(dat.1$Treatment, levels = c("Vehicle", "Amp", "Metro", "Amp + Metro"))
levels(dat.1$Treatment)

## [1] "Vehicle"     "Amp"         "Metro"       "Amp + Metro"

p.boxplot.phylum.2 <- ggpaired(dat.1, x = "DaysTreatment", y = "Abundance", outlier.shape = NA, id = "CageDuringTreatment", line.size = 0.5, line.color = "lightgray", point.size = 0.3, width = 0.4) +
  #geom_point(size = 1) +
  #geom_jitter(width = 0.2, alpha = 0.7) +
  ylab("Relative Abundance") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  stat_compare_means(label = "p.signif", method = "t.test", ref.group = "D0", hide.ns = TRUE, label.y.npc = 0.95) +
  facet_grid(Treatment~Phylum) +
  theme(axis.text.x = element_text(size = 6)) +
  theme(axis.text.y = element_text(size = 6)) +
  theme(axis.title.y = element_text(size = 6)) +
  theme(axis.title.x = element_blank()) +
  theme(strip.text = element_text(size = 6)) +
  ylim(0,1)
p.boxplot.phylum.2

## Phylum level general additive model (GAM) plots

# Define color scheme
my.cols <- brewer.pal(n = 8, "Dark2")
my.cols[3] <- "#08519C"

# Phyla plots with GAM smoother 
p.gam.phylum <- ggplot(dat.1, aes(x = Day, y = Abundance, color = Phylum, group = Phylum)) +
  stat_smooth(method = "gam", formula = y ~ s(x, bs = "cr", k = 7)) +
  facet_grid(~Treatment) +
  ylab("Relative Abundance") +
  geom_point(size = 1.25, alpha = 0.4) +
  theme(axis.text.x = element_text(size = 10)) +
  theme(axis.text.y = element_text(size = 10)) +
  theme(axis.title.y = element_text(size = 10)) +
  theme(axis.title.x = element_blank()) +
  theme(strip.text = element_text(size = 10)) +
  scale_y_continuous(breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2)) +
  scale_x_continuous(breaks = c(0, 3, 7, 13, 16, 18, 20)) +
  scale_color_manual(values = my.cols) +
  theme(strip.background = element_blank()) +
  theme(strip.text.x = element_blank()) +
  theme(axis.title.y = element_blank()) +
  guides(color=guide_legend(override.aes=list(fill=NA))) +
  theme(legend.title = element_blank())
p.gam.phylum

# Testing plots - Close up view of how each phyla behaves
dat.bac <- subset(dat, Phylum == "Bacteroidetes")
dat.firm <- subset(dat, Phylum == "Firmicutes")
dat.teneri <- subset(dat, Phylum == "Tenericutes")
dat.proteo <- subset(dat, Phylum == "Proteobacteria")

Phylum level general additive model (GAM) analysis

# Convert data frame so that each animal (CageDuringTreatment) has a unique phyla and abundnace (long form transformation with melt)
dat.1.melt <- melt(dat.1, id.vars = c("CageDuringTreatment", "Day", "Phylum", "Treatment"), measure.vars = c("Abundance"))

# Recast the data by indiviudal phyla
dat.1.cast <- dcast(dat.1.melt, CageDuringTreatment + Treatment + Day ~ Phylum)

# Subset melt/casted frames by Treatment for within treatment testing
dat.1.cast.vehicle <- filter(dat.1.cast, Treatment %in% "Vehicle")
dat.1.cast.amp <- filter(dat.1.cast, Treatment %in% "Amp")
dat.1.cast.metro <- filter(dat.1.cast, Treatment %in% "Metro")
dat.1.cast.ampmetro <- filter(dat.1.cast, Treatment %in% "Amp + Metro")

## General Additive Model (GAM) testing for the contribution of the relative abundance of each Phyla to the change in abundnace over time
# Vehicle
mod_gam.vehicle.1 <- gam(Day ~ s(Bacteroidetes, bs = "cr", k = 7) + s(Firmicutes, bs = "cr", k = 7) + s(Proteobacteria, bs = "cr", k = 7) + s(Tenericutes, bs = "cr", k = 7), data = dat.1.cast.vehicle)
summary(mod_gam.vehicle.1)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Day ~ s(Bacteroidetes, bs = "cr", k = 7) + s(Firmicutes, bs = "cr", 
##     k = 7) + s(Proteobacteria, bs = "cr", k = 7) + s(Tenericutes, 
##     bs = "cr", k = 7)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   11.000      0.742   14.82   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                      edf Ref.df     F p-value  
## s(Bacteroidetes)  0.9999  1.000 3.823  0.0550 .
## s(Firmicutes)     1.0000  1.000 4.377  0.0404 *
## s(Proteobacteria) 2.0598  2.473 2.946  0.0457 *
## s(Tenericutes)    2.5961  3.102 1.925  0.1226  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.261   Deviance explained = 33.3%
## GCV = 43.271  Scale est. = 38.539    n = 70

ggpaired(subset(dat.proteo, Treatment == "Vehicle"), x = "DaysTreatment", y = "Abundance", outlier.shape = NA, id = "CageDuringTreatment", line.size = 0.5, line.color = "lightgray", point.size = 0.3, width = 0.4) +
  stat_compare_means(label = "p.signif", method = "t.test", ref.group = "D0", hide.ns = TRUE, label.y.npc = 1) +
  ggtitle("Proteobacteria in Vehicle Control Mice")

# Metro
mod_gam.metro.1 <- gam(Day ~ s(Bacteroidetes, bs = "cr", k = 7) + s(Firmicutes, bs = "cr", k = 7) + s(Proteobacteria, bs = "cr", k = 7) + s(Tenericutes, bs = "cr", k = 7), data = dat.1.cast.metro)
summary(mod_gam.metro.1)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Day ~ s(Bacteroidetes, bs = "cr", k = 7) + s(Firmicutes, bs = "cr", 
##     k = 7) + s(Proteobacteria, bs = "cr", k = 7) + s(Tenericutes, 
##     bs = "cr", k = 7)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  11.0000     0.6067   18.13   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                     edf Ref.df     F p-value    
## s(Bacteroidetes)  5.467  5.796 1.710  0.1949    
## s(Firmicutes)     4.314  5.084 1.582  0.1809    
## s(Proteobacteria) 4.694  5.250 2.766  0.0205 *  
## s(Tenericutes)    5.809  5.969 8.241 5.4e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.506   Deviance explained = 65.1%
## GCV = 37.023  Scale est. = 25.766    n = 70

ggpaired(subset(dat.teneri, Treatment == "Metro"), x = "DaysTreatment", y = "Abundance", outlier.shape = NA, id = "CageDuringTreatment", line.size = 0.5, line.color = "lightgray", point.size = 0.3, width = 0.4) +
  stat_compare_means(label = "p.signif", method = "t.test", ref.group = "D0", hide.ns = TRUE, label.y.npc = 1) +
  ggtitle("Tenericutes in Metro Treated Mice")

ggpaired(subset(dat.proteo, Treatment == "Metro"), x = "DaysTreatment", y = "Abundance", outlier.shape = NA, id = "CageDuringTreatment", line.size = 0.5, line.color = "lightgray", point.size = 0.3, width = 0.4) +
  stat_compare_means(label = "p.signif", method = "t.test", ref.group = "D0", hide.ns = TRUE, label.y.npc = 1) +
  ggtitle("Proteobacteria in Metro Treated Mice")

# Amp
mod_gam.amp.1 <- gam(Day ~ s(Bacteroidetes, bs = "cr", k = 7) + s(Firmicutes, bs = "cr", k = 7) + s(Proteobacteria, bs = "cr", k = 7) + s(Tenericutes, bs = "cr", k = 7), data = dat.1.cast.amp)
summary(mod_gam.amp.1)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Day ~ s(Bacteroidetes, bs = "cr", k = 7) + s(Firmicutes, bs = "cr", 
##     k = 7) + s(Proteobacteria, bs = "cr", k = 7) + s(Tenericutes, 
##     bs = "cr", k = 7)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  11.0000     0.5499      20   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                     edf Ref.df     F  p-value    
## s(Bacteroidetes)  1.000  1.000 0.193 0.661842    
## s(Firmicutes)     5.023  5.596 3.057 0.016504 *  
## s(Proteobacteria) 5.264  5.632 3.437 0.009471 ** 
## s(Tenericutes)    4.725  5.007 4.598 0.000863 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.388   Deviance explained = 48.3%
## GCV =  37.89  Scale est. = 31.751    n = 105

# Amp + Metro
mod_gam.ampmetro <- gam(Day ~ s(Bacteroidetes, bs = "cr", k = 7) + s(Firmicutes, bs = "cr", k = 7) + s(Proteobacteria, bs = "cr", k = 7) + s(Tenericutes, bs = "cr", k = 7), data = dat.1.cast.ampmetro)
summary(mod_gam.ampmetro)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Day ~ s(Bacteroidetes, bs = "cr", k = 7) + s(Firmicutes, bs = "cr", 
##     k = 7) + s(Proteobacteria, bs = "cr", k = 7) + s(Tenericutes, 
##     bs = "cr", k = 7)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  10.3438     0.5071    20.4   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                      edf Ref.df       F  p-value    
## s(Bacteroidetes)  1.0544 1.2899  28.445 0.000270 ***
## s(Firmicutes)     0.7196 0.7678   0.065 0.824564    
## s(Proteobacteria) 0.4675 0.4717 140.607 3.39e-12 ***
## s(Tenericutes)    2.3332 2.6542   7.372 0.000382 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 23/25
## R-sq.(adj) =  0.682   Deviance explained = 70.5%
## GCV = 18.031  Scale est. = 16.46     n = 64

Alpha diversity summary information generation

# Diversity
diversity <- global(ps1)
head(diversity)

##                    richness_0 richness_20 richness_50 richness_80
## 101.Thackray.D0.E1        121         121         121         121
## 102.Thackray.D0.E2        166         166         166         166
## 103.Thackray.D0.E3         93          93          93          93
## 104.Thackray.D0.E4        151         151         151         151
## 105.Thackray.D0.E5        173         173         173         173
## 106.Thackray.D0.F1        102         102         102         102
##                    diversities_inverse_simpson diversities_gini_simpson
## 101.Thackray.D0.E1                   12.170885                0.9178367
## 102.Thackray.D0.E2                   25.401667                0.9606325
## 103.Thackray.D0.E3                   10.355095                0.9034292
## 104.Thackray.D0.E4                   24.015585                0.9583604
## 105.Thackray.D0.E5                   30.875326                0.9676117
## 106.Thackray.D0.F1                    9.038626                0.8893637
##                    diversities_shannon diversities_fisher
## 101.Thackray.D0.E1            3.081963           16.60700
## 102.Thackray.D0.E2            3.983350           23.16563
## 103.Thackray.D0.E3            2.923112           12.47155
## 104.Thackray.D0.E4            3.802921           21.14357
## 105.Thackray.D0.E5            4.068889           24.48143
## 106.Thackray.D0.F1            2.836865           13.94334
##                    diversities_coverage evenness_camargo evenness_pielou
## 101.Thackray.D0.E1                    5       0.01410119       0.6426391
## 102.Thackray.D0.E2                   10       0.03354858       0.7792174
## 103.Thackray.D0.E3                    4       0.01148292       0.6449084
## 104.Thackray.D0.E4                    8       0.02752578       0.7579648
## 105.Thackray.D0.E5                   11       0.03541365       0.7895708
## 106.Thackray.D0.F1                    3       0.01154458       0.6133797
##                    evenness_simpson evenness_evar evenness_bulla
## 101.Thackray.D0.E1      0.008929483     0.2255497     0.06979504
## 102.Thackray.D0.E2      0.018636586     0.2804325     0.10587282
## 103.Thackray.D0.E3      0.007597282     0.2026502     0.05605963
## 104.Thackray.D0.E4      0.017619651     0.2691205     0.09552794
## 105.Thackray.D0.E5      0.022652477     0.2584945     0.10681799
## 106.Thackray.D0.F1      0.006631420     0.2367310     0.06282215
##                    dominance_dbp dominance_dmn dominance_absolute
## 101.Thackray.D0.E1    0.16531143     0.2888926               4005
## 102.Thackray.D0.E2    0.13306855     0.2103331               3987
## 103.Thackray.D0.E3    0.20086163     0.3593830               4336
## 104.Thackray.D0.E4    0.09593916     0.1840488               2561
## 105.Thackray.D0.E5    0.10275191     0.1638590               2946
## 106.Thackray.D0.F1    0.19345906     0.3746479               4052
##                    dominance_relative dominance_simpson
## 101.Thackray.D0.E1         0.16531143        0.08216330
## 102.Thackray.D0.E2         0.13306855        0.03936749
## 103.Thackray.D0.E3         0.20086163        0.09657082
## 104.Thackray.D0.E4         0.09593916        0.04163963
## 105.Thackray.D0.E5         0.10275191        0.03238832
## 106.Thackray.D0.F1         0.19345906        0.11063629
##                    dominance_core_abundance dominance_gini
## 101.Thackray.D0.E1                0.7752920      0.9858988
## 102.Thackray.D0.E2                0.5431547      0.9664514
## 103.Thackray.D0.E3                0.7725483      0.9885171
## 104.Thackray.D0.E4                0.6345995      0.9724742
## 105.Thackray.D0.E5                0.4216804      0.9645863
## 106.Thackray.D0.F1                0.8083075      0.9884554
##                    rarity_log_modulo_skewness rarity_low_abundance
## 101.Thackray.D0.E1                   2.059941           0.05733273
## 102.Thackray.D0.E2                   2.058710           0.07970095
## 103.Thackray.D0.E3                   2.059384           0.04660212
## 104.Thackray.D0.E4                   2.059367           0.06799281
## 105.Thackray.D0.E5                   2.060132           0.06867567
## 106.Thackray.D0.F1                   2.060471           0.05743614
##                    rarity_noncore_abundance rarity_rare_abundance
## 101.Thackray.D0.E1               0.08527676            0.03178272
## 102.Thackray.D0.E2               0.27635004            0.09892531
## 103.Thackray.D0.E3               0.10779636            0.04377635
## 104.Thackray.D0.E4               0.19850903            0.05952649
## 105.Thackray.D0.E5               0.28607304            0.11677305
## 106.Thackray.D0.F1               0.10188589            0.06378611

ps1.rich <- merge(sd, diversity, by ="row.names") # merge sd.1 by row names

# Add divergence measurements
ps1.rich$divergence <- divergence(ps1)

Alpha diversity analysis - Paired Boxplots

p.rich.treatment <- ggpaired(ps1.rich,  x = "DaysTreatment", y = "richness_0", outlier.shape = NA, id = "CageDuringTreatment", line.size = 0.5, line.color = "gray") +
  geom_jitter(width = 0.2) +
  facet_grid(~Treatment) +
  ylab("Richness") +
  theme(axis.title.x = element_blank()) +
  stat_compare_means(label = "p.signif", method = "t.test", ref.group = "D0", hide.ns = TRUE)

p.sd.treatment <- ggpaired(ps1.rich, x = "DaysTreatment", y = "diversities_shannon", outlier.shape = NA, id = "CageDuringTreatment", line.size = 0.5, line.color = "gray") +
  geom_boxplot(position = position_dodge()) +
  geom_jitter(width = 0.2) +
  facet_grid(~Treatment) +
  ylab("Shannon diversity") +
  theme(axis.title.x = element_blank()) +
  stat_compare_means(label = "p.signif", method = "t.test", ref.group = "D0", hide.ns = TRUE)

grid.arrange(p.rich.treatment, p.sd.treatment, nrow = 2)

## Alpha diversity general additive model (GAM) analysis.

ps1.rich.melt <- melt(ps1.rich, id.vars = c("Treatment", "Day", "DaysTreatment"), measure.vars = c("richness_0"))
ps1.sd.melt <- melt(ps1.rich, id.vars = c("Treatment", "Day", "DaysTreatment"), measure.vars = c("diversities_shannon"))

# Richness
p.rich.gam.treat <- ggplot(ps1.rich.melt, aes(x = Day, y = value, color = Treatment, group = Treatment)) +
  stat_smooth(method = "gam", formula = y ~ s(x, bs = "cr", k = 7)) +
  ylab("Richness") +
  geom_point(size = 1.25, alpha = 0.5) +
  theme(axis.text.x = element_text(size = 10)) +
  theme(axis.text.y = element_text(size = 10)) +
  theme(axis.title.y = element_blank()) +
  theme(axis.title.x = element_blank()) +
  theme(strip.text = element_text(size = 10)) +
  scale_color_manual(values = c("black", "chocolate", "green", "purple")) +
  theme(legend.position = "NULL") +
  scale_y_continuous(limits = c(0,250), breaks = c(0,50, 100, 150, 200, 250)) +
  scale_x_continuous(breaks = c(0,3,7,13,16,18,20))

# Shannon diversity
p.sd.gam.treat <- ggplot(ps1.sd.melt, aes(x = Day, y = value, color = Treatment, group = Treatment)) +
  stat_smooth(method = "gam", formula = y ~ s(x, bs = "cr", k = 7)) +
  ylab("Shannon Diversity") +
  geom_point(size = 1.25, alpha = 0.5) +
  theme(axis.text.x = element_text(size = 10)) +
  theme(axis.text.y = element_text(size = 10)) +
  theme(axis.title.y = element_blank()) +
  theme(axis.title.x = element_blank()) +
  theme(strip.text = element_text(size = 10)) +
  scale_color_manual(values = c("black", "green", "chocolate", "purple"), labels = c("Vehicle", "A", "M", "AM")) +
  guides(color=guide_legend(override.aes=list(fill=NA))) +
  theme(legend.position = "right") +
  theme(legend.title = element_blank()) +
  scale_y_continuous(limits = c(0,5), breaks = c(0,1,2,3,4,5)) +
  scale_x_continuous(breaks = c(0,3,7,13,16,18,20))
  #theme(legend.title = element_blank())

grid.arrange(p.rich.gam.treat, p.sd.gam.treat, ncol = 2)

## Ordination

#Ordination Analysis
#Beta Diversity has same trend of timepoints with longtail and bimodal read counts having larger elipses
ord.pcoa.bray <- ordinate(ps1, method = "PCoA", distance = "bray")
ord.pcoa.uni <- ordinate(ps1, method = "PCoA", distance = "unifrac")
ord.pcoa.wuni <- ordinate(ps1, method = "PCoA", distance = "wunifrac")

Beta diversity ordination plots ~ SurvivalStatus

## Ordination plots all samples
# Bray
p.pcoa.bray <- plot_ordination(ps1, ord.pcoa.bray, color = "SurvivalStatus", axes = c(1,2)) +
  geom_point(size = 2) +
  ggtitle("PCoA of UniFrac Distances") +
  facet_grid(Treatment~DaysTreatment)
  #stat_ellipse(type = "norm", geom = "polygon", alpha = 1/10, aes(fill = SurvivalStatus))
p.pcoa.bray

# Unifrac
p.pcoa.uni <- plot_ordination(ps1, ord.pcoa.uni, color = "SurvivalStatus", axes = c(1,2)) +
  geom_point(size = 2) +
  ggtitle("PCoA of UniFrac Distances") +
  facet_grid(Treatment~DaysTreatment)
  #stat_ellipse(type = "norm", geom = "polygon", alpha = 1/10, aes(fill = SurvivalStatus))
p.pcoa.uni

# Weighted Unifrac
p.pcoa.wuni <- plot_ordination(ps1, ord.pcoa.wuni, color = "SurvivalStatus", axes = c(1,2)) +
  geom_point(size = 2) +
  ggtitle("PCoA of wUniFrac Distances") +
  facet_grid(Treatment~DaysTreatment) +
  stat_ellipse(type = "norm", geom = "polygon", alpha = 1/10, aes(fill = SurvivalStatus))
p.pcoa.wuni

p.pcoa.uni.treat <- plot_ordination(ps1, ord.pcoa.uni, color = "Treatment", shape = "SurvivalStatus") +
  geom_point(size = 3) +
  # ggtitle("PCoA of UniFrac Distances") +
  facet_grid(~DaysTreatment) +
  scale_color_manual(values = c("black", "green", "chocolate", "purple"), labels = c("Vehicle", "A", "M", "AM")) +
  theme(axis.title.y = element_blank()) +
  theme(axis.title.x = element_blank()) +
  theme(legend.position = "bottom") +
  theme(legend.text = element_text(size = 10)) +
  theme(legend.title = element_blank()) +
  theme(strip.background = element_blank()) +
  theme(strip.text.x = element_blank())
p.pcoa.uni.treat

##Group significance testing with ADONIS

# Set a random seed so that exact results can be reproduced
set.seed(10000)

# Function to run adonis test on a physeq object and a variable from metadata 
doadonis <- function(physeq, category) {
  bdist <- phyloseq::distance(physeq, "unifrac")
  col <- as(sample_data(physeq), "data.frame")[ ,category]
  
  # Adonis test
  adonis.bdist <- adonis(bdist ~ col)
  print("Adonis results:")
  print(adonis.bdist)
  
  # Homogeneity of dispersion test
  betatax = betadisper(bdist,col)
  p = permutest(betatax)
  print("Betadisper results:")
  print(p$tab)
}

doadonis(ps1, "Treatment")

## [1] "Adonis results:"
## 
## Call:
## adonis(formula = bdist ~ col) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##            Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## col         3    17.812  5.9375  39.413 0.27937  0.001 ***
## Residuals 305    45.947  0.1506         0.72063           
## Total     308    63.760                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Betadisper results:"
##            Df   Sum Sq     Mean Sq        F N.Perm Pr(>F)
## Groups      3 4.419388 1.473129306 262.2619    999  0.001
## Residuals 305 1.713190 0.005617016       NA     NA     NA

library(pairwiseAdonis)

ps1.d0 <- subset_samples(ps1, DaysTreatment == "D0")
ps1_otu_table.d0 <- as.data.frame(otu_table(ps1.d0))
sd.df.d0 <- as.data.frame(sample_data(ps1.d0))
ps1_otu_table.d0$Treatment <- sd.df.d0$Treatment

ps1.d3 <- subset_samples(ps1, DaysTreatment == "D3")
ps1_otu_table.d3 <- as.data.frame(otu_table(ps1.d3))
sd.df.d3 <- as.data.frame(sample_data(ps1.d3))
ps1_otu_table.d3$Treatment <- sd.df.d3$Treatment

ps1.d7 <- subset_samples(ps1, DaysTreatment == "D7")
ps1_otu_table.d7 <- as.data.frame(otu_table(ps1.d7))
sd.df.d7 <- as.data.frame(sample_data(ps1.d7))
ps1_otu_table.d7$Treatment <- sd.df.d7$Treatment

ps1.d13 <- subset_samples(ps1, DaysTreatment == "D13")
ps1_otu_table.d13 <- as.data.frame(otu_table(ps1.d13))
sd.df.d13 <- as.data.frame(sample_data(ps1.d13))
ps1_otu_table.d13$Treatment <- sd.df.d13$Treatment

ps1.d16 <- subset_samples(ps1, DaysTreatment == "D16")
ps1_otu_table.d16 <- as.data.frame(otu_table(ps1.d16))
sd.df.d16 <- as.data.frame(sample_data(ps1.d16))
ps1_otu_table.d16$Treatment <- sd.df.d16$Treatment

ps1.d20 <- subset_samples(ps1, DaysTreatment == "D20")
ps1_otu_table.d20 <- as.data.frame(otu_table(ps1.d20))
sd.df.d20 <- as.data.frame(sample_data(ps1.d20))
ps1_otu_table.d20$Treatment <- sd.df.d20$Treatment

pairwise.adonis(ps1_otu_table.d0[,1:1363], ps1_otu_table.d0$Treatment)

## [1] "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1"

##                    pairs  F.Model         R2 p.value p.adjusted sig
## 1           Metro vs Amp 2.541514 0.09950522   0.028      0.168    
## 2   Metro vs Amp + Metro 3.012322 0.14335980   0.008      0.048   .
## 3       Metro vs Vehicle 3.008221 0.14319254   0.020      0.120    
## 4     Amp vs Amp + Metro 0.952873 0.03978116   0.465      1.000    
## 5         Amp vs Vehicle 2.675326 0.10419834   0.010      0.060    
## 6 Amp + Metro vs Vehicle 3.143583 0.14867786   0.001      0.006   *

pairwise.adonis(ps1_otu_table.d3[,1:1363], ps1_otu_table.d3$Treatment)

## [1] "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1"

##                    pairs   F.Model        R2 p.value p.adjusted sig
## 1       Vehicle vs Metro  8.531852 0.3215702   0.001      0.006   *
## 2         Vehicle vs Amp 17.429204 0.4311043   0.001      0.006   *
## 3 Vehicle vs Amp + Metro  8.652032 0.3246294   0.001      0.006   *
## 4           Metro vs Amp 20.688576 0.4735466   0.001      0.006   *
## 5   Metro vs Amp + Metro 11.193323 0.3834207   0.001      0.006   *
## 6     Amp vs Amp + Metro  3.885616 0.1445240   0.001      0.006   *

pairwise.adonis(ps1_otu_table.d7[,1:1363], ps1_otu_table.d7$Treatment)

## [1] "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1"

##                    pairs   F.Model        R2 p.value p.adjusted sig
## 1       Vehicle vs Metro  5.590241 0.2369726   0.001      0.006   *
## 2         Vehicle vs Amp 18.023438 0.4393449   0.001      0.006   *
## 3 Vehicle vs Amp + Metro 11.673062 0.3933892   0.001      0.006   *
## 4           Metro vs Amp 27.042017 0.5403862   0.001      0.006   *
## 5   Metro vs Amp + Metro 16.419568 0.4770417   0.001      0.006   *
## 6     Amp vs Amp + Metro  4.396955 0.1604906   0.001      0.006   *

pairwise.adonis(ps1_otu_table.d13[,1:1363], ps1_otu_table.d13$Treatment)

## [1] "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1"

##                    pairs   F.Model        R2 p.value p.adjusted sig
## 1       Vehicle vs Metro  3.723724 0.1714128   0.002      0.012   .
## 2         Vehicle vs Amp 19.990937 0.4650035   0.001      0.006   *
## 3 Vehicle vs Amp + Metro 18.657895 0.5089734   0.001      0.006   *
## 4           Metro vs Amp 20.031006 0.4655017   0.001      0.006   *
## 5   Metro vs Amp + Metro 19.103878 0.5148755   0.001      0.006   *
## 6     Amp vs Amp + Metro  9.760540 0.2979359   0.001      0.006   *

pairwise.adonis(ps1_otu_table.d16[,1:1363], ps1_otu_table.d16$Treatment)

## [1] "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1"

##                    pairs   F.Model        R2 p.value p.adjusted sig
## 1       Vehicle vs Metro  2.568894 0.1248922   0.014      0.084    
## 2         Vehicle vs Amp 20.190996 0.4674816   0.001      0.006   *
## 3 Vehicle vs Amp + Metro 32.948284 0.6731244   0.001      0.006   *
## 4           Metro vs Amp 16.937747 0.4241037   0.001      0.006   *
## 5   Metro vs Amp + Metro 24.788139 0.6077291   0.001      0.006   *
## 6     Amp vs Amp + Metro 27.745340 0.5691896   0.001      0.006   *

pairwise.adonis(ps1_otu_table.d20[,1:1363], ps1_otu_table.d20$Treatment)

## [1] "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1"

##                    pairs   F.Model        R2 p.value p.adjusted sig
## 1       Vehicle vs Metro   6.82161 0.2748254   0.001      0.006   *
## 2         Vehicle vs Amp  16.74809 0.4213559   0.001      0.006   *
## 3 Vehicle vs Amp + Metro  47.10367 0.7464490   0.001      0.006   *
## 4           Metro vs Amp  21.72066 0.4856963   0.001      0.006   *
## 5   Metro vs Amp + Metro 109.84593 0.8728604   0.001      0.006   *
## 6     Amp vs Amp + Metro  19.88933 0.4864185   0.001      0.006   *

# Amp ~ Day
ps1.amp.d0 <- subset_samples(ps1.amp, DaysTreatment == "D0")
ps1.amp.d3 <- subset_samples(ps1.amp, DaysTreatment == "D3")
ps1.amp.d7 <- subset_samples(ps1.amp, DaysTreatment == "D7")
ps1.amp.d13 <- subset_samples(ps1.amp, DaysTreatment == "D13")
ps1.amp.d16 <- subset_samples(ps1.amp, DaysTreatment == "D16")
ps1.amp.d20 <- subset_samples(ps1.amp, DaysTreatment == "D20")

doadonis(ps1.amp.d0, "SurvivalStatus") # *

## [1] "Adonis results:"
## 
## Call:
## adonis(formula = bdist ~ col) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs  MeanSqs F.Model      R2 Pr(>F)
## col        1   0.05114 0.051141  1.0061 0.07183  0.418
## Residuals 13   0.66079 0.050830         0.92817       
## Total     14   0.71193                  1.00000       
## [1] "Betadisper results:"
##           Df    Sum Sq     Mean Sq        F N.Perm Pr(>F)
## Groups     1 0.0040013 0.004001300 2.693166    999  0.117
## Residuals 13 0.0193144 0.001485723       NA     NA     NA

doadonis(ps1.amp.d3, "SurvivalStatus") # n.s.

## [1] "Adonis results:"
## 
## Call:
## adonis(formula = bdist ~ col) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
## col        1   0.15533 0.15533 0.83339 0.06024  0.647
## Residuals 13   2.42299 0.18638         0.93976       
## Total     14   2.57832                 1.00000       
## [1] "Betadisper results:"
##           Df      Sum Sq     Mean Sq          F N.Perm Pr(>F)
## Groups     1 0.000335951 0.000335951 0.01707187    999  0.917
## Residuals 13 0.255822190 0.019678630         NA     NA     NA

doadonis(ps1.amp.d7, "SurvivalStatus") # n.s.

## [1] "Adonis results:"
## 
## Call:
## adonis(formula = bdist ~ col) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
## col        1   0.13571 0.13572  1.1761 0.08297   0.25
## Residuals 13   1.50008 0.11539         0.91703       
## Total     14   1.63579                 1.00000       
## [1] "Betadisper results:"
##           Df      Sum Sq     Mean Sq        F N.Perm Pr(>F)
## Groups     1 0.007527834 0.007527834 2.554168    999  0.124
## Residuals 13 0.038314564 0.002947274       NA     NA     NA

doadonis(ps1.amp.d13, "SurvivalStatus") # n.s. 

## [1] "Adonis results:"
## 
## Call:
## adonis(formula = bdist ~ col) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
## col        1   0.21715 0.21715  1.4577 0.10083  0.136
## Residuals 13   1.93658 0.14897         0.89917       
## Total     14   2.15374                 1.00000       
## [1] "Betadisper results:"
##           Df     Sum Sq     Mean Sq        F N.Perm Pr(>F)
## Groups     1 0.02576481 0.025764813 4.057142    999  0.075
## Residuals 13 0.08255628 0.006350483       NA     NA     NA

doadonis(ps1.amp.d16, "SurvivalStatus") # n.s.

## [1] "Adonis results:"
## 
## Call:
## adonis(formula = bdist ~ col) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
## col        1   0.20463 0.20463  1.0978 0.07787  0.328
## Residuals 13   2.42313 0.18639         0.92213       
## Total     14   2.62776                 1.00000       
## [1] "Betadisper results:"
##           Df     Sum Sq     Mean Sq        F N.Perm Pr(>F)
## Groups     1 0.01037965 0.010379646 1.535046    999  0.213
## Residuals 13 0.08790315 0.006761781       NA     NA     NA

doadonis(ps1.amp.d20, "SurvivalStatus") # n.s.

## [1] "Adonis results:"
## 
## Call:
## adonis(formula = bdist ~ col) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
## col        1    0.1373 0.13734 0.58365 0.04297  0.887
## Residuals 13    3.0591 0.23531         0.95703       
## Total     14    3.1964                 1.00000       
## [1] "Betadisper results:"
##           Df      Sum Sq     Mean Sq         F N.Perm Pr(>F)
## Groups     1 0.001757289 0.001757289 0.1385605    999  0.695
## Residuals 13 0.164872085 0.012682468        NA     NA     NA

# Treatment ~ Day
ps1.d0 <- subset_samples(ps1, DaysTreatment == "D0")
ps1.d3 <- subset_samples(ps1, DaysTreatment == "D3")
ps1.d7 <- subset_samples(ps1, DaysTreatment == "D7")
ps1.d13 <- subset_samples(ps1, DaysTreatment == "D13")
ps1.d16 <- subset_samples(ps1, DaysTreatment == "D16")
ps1.d20 <- subset_samples(ps1, DaysTreatment == "D20")

doadonis(ps1.d0, "Treatment") # *

## [1] "Adonis results:"
## 
## Call:
## adonis(formula = bdist ~ col) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs  MeanSqs F.Model      R2 Pr(>F)    
## col        3   0.36198 0.120660  2.2241 0.13996  0.001 ***
## Residuals 41   2.22432 0.054252         0.86004           
## Total     44   2.58630                  1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Betadisper results:"
##           Df      Sum Sq     Mean Sq       F N.Perm Pr(>F)
## Groups     3 0.006101085 0.002033695 1.53105    999  0.231
## Residuals 41 0.054460352 0.001328301      NA     NA     NA

doadonis(ps1.d3, "Treatment") # **

## [1] "Adonis results:"
## 
## Call:
## adonis(formula = bdist ~ col) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## col        3    3.5719 1.19063  9.8536 0.41894  0.001 ***
## Residuals 41    4.9541 0.12083         0.58106           
## Total     44    8.5260                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Betadisper results:"
##           Df    Sum Sq     Mean Sq        F N.Perm Pr(>F)
## Groups     3 0.3083364 0.102778794 11.34675    999  0.001
## Residuals 41 0.3713778 0.009057995       NA     NA     NA

doadonis(ps1.d7, "Treatment") # ***

## [1] "Adonis results:"
## 
## Call:
## adonis(formula = bdist ~ col) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## col        3    3.1929 1.06429  11.477 0.45645  0.001 ***
## Residuals 41    3.8021 0.09273         0.54355           
## Total     44    6.9949                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Betadisper results:"
##           Df    Sum Sq     Mean Sq        F N.Perm Pr(>F)
## Groups     3 0.3818544 0.127284804 33.01172    999  0.001
## Residuals 41 0.1580856 0.003855746       NA     NA     NA

doadonis(ps1.d13, "Treatment") # ***

## [1] "Adonis results:"
## 
## Call:
## adonis(formula = bdist ~ col) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## col        3    4.1878  1.3959  11.498 0.45692  0.001 ***
## Residuals 41    4.9775  0.1214         0.54308           
## Total     44    9.1653                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Betadisper results:"
##           Df    Sum Sq     Mean Sq        F N.Perm Pr(>F)
## Groups     3 0.3405377 0.113512572 14.31893    999  0.001
## Residuals 41 0.3250254 0.007927448       NA     NA     NA

doadonis(ps1.d16, "Treatment") # ***

## [1] "Adonis results:"
## 
## Call:
## adonis(formula = bdist ~ col) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## col        3    5.3900 1.79667  14.644 0.52974  0.001 ***
## Residuals 39    4.7849 0.12269         0.47026           
## Total     42   10.1749                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Betadisper results:"
##           Df    Sum Sq     Mean Sq        F N.Perm Pr(>F)
## Groups     3 0.3688410 0.122946985 17.11128    999  0.001
## Residuals 39 0.2802205 0.007185142       NA     NA     NA

doadonis(ps1.d20, "Treatment") # ***

## [1] "Adonis results:"
## 
## Call:
## adonis(formula = bdist ~ col) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## col        3    5.2036 1.73452  15.837 0.54919  0.001 ***
## Residuals 39    4.2714 0.10952         0.45081           
## Total     42    9.4750                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Betadisper results:"
##           Df    Sum Sq     Mean Sq        F N.Perm Pr(>F)
## Groups     3 0.8190444 0.273014799 41.90247    999  0.001
## Residuals 39 0.2541038 0.006515482       NA     NA     NA

NMDS time-series analysis

R session info

sessionInfo()

## R version 3.4.3 (2017-11-30)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.2
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] bindrcpp_0.2         pairwiseAdonis_0.0.1 cluster_2.0.6       
##  [4] mgcv_1.8-22          nlme_3.1-131         data.table_1.10.4-3 
##  [7] ggpubr_0.1.6         magrittr_1.5         microbiome_1.0.0    
## [10] plotly_4.7.1         knitr_1.18           gridExtra_2.3       
## [13] vegan_2.4-5          lattice_0.20-35      permute_0.9-4       
## [16] RColorBrewer_1.1-2   phyloseq_1.23.1      plyr_1.8.4          
## [19] reshape2_1.4.3       forcats_0.2.0        stringr_1.2.0       
## [22] dplyr_0.7.4          purrr_0.2.4          readr_1.1.1         
## [25] tidyr_0.7.2          tibble_1.4.1         ggplot2_2.2.1.9000  
## [28] tidyverse_1.2.1     
## 
## loaded via a namespace (and not attached):
##  [1] lubridate_1.7.1     httr_1.3.1          rprojroot_1.3-2    
##  [4] tools_3.4.3         backports_1.1.2     R6_2.2.2           
##  [7] lazyeval_0.2.1      BiocGenerics_0.24.0 colorspace_1.3-2   
## [10] ade4_1.7-10         withr_2.1.1         tidyselect_0.2.3   
## [13] mnormt_1.5-5        compiler_3.4.3      cli_1.0.0          
## [16] rvest_0.3.2         Biobase_2.38.0      xml2_1.1.1         
## [19] labeling_0.3        scales_0.5.0.9000   psych_1.7.8        
## [22] digest_0.6.14       foreign_0.8-69      rmarkdown_1.8      
## [25] XVector_0.18.0      pkgconfig_2.0.1     htmltools_0.3.6    
## [28] htmlwidgets_0.9     rlang_0.1.6         readxl_1.0.0       
## [31] rstudioapi_0.7      shiny_1.0.5         bindr_0.1          
## [34] jsonlite_1.5        crosstalk_1.0.0     biomformat_1.6.0   
## [37] Matrix_1.2-12       Rcpp_0.12.14        munsell_0.4.3      
## [40] S4Vectors_0.16.0    ape_5.0             stringi_1.1.6      
## [43] yaml_2.1.16         MASS_7.3-48         zlibbioc_1.24.0    
## [46] rhdf5_2.22.0        grid_3.4.3          parallel_3.4.3     
## [49] crayon_1.3.4        Biostrings_2.46.0   haven_1.1.0        
## [52] splines_3.4.3       multtest_2.34.0     hms_0.4.0          
## [55] pillar_1.0.1        igraph_1.1.2        codetools_0.2-15   
## [58] stats4_3.4.3        glue_1.2.0          evaluate_0.10.1    
## [61] modelr_0.1.1        httpuv_1.3.5        foreach_1.4.4      
## [64] cellranger_1.1.0    gtable_0.2.0        assertthat_0.2.0   
## [67] mime_0.5            xtable_1.8-2        broom_0.4.3        
## [70] survival_2.41-3     viridisLite_0.2.0   iterators_1.0.9    
## [73] IRanges_2.12.0